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ABSTRACT 


iQyf}p\Qxi-yc-  tr 


Reactive  systems  involving  Arrhenius  kinetics  often  exhibit  multiple 
steady  states.  A  typical  response  is  an  S-curve,  whose  turnaround  points 
correspond  to  ignition  or  extinction.  This  paper  describes  the  dynamics  of 
transition  from  the  extinguished  to  the  ignited  state  as  the  reaction-rate 
parameter  is  slowly  varied  through  the  critical  value.  Both  lumped  and 
spatially  distributed  models  are  studied.  The  asymptotic  analysis  is  based 
on  the  largeness  of  two  parameters:  one  characterizing  the  activation 
energy  and  the  other  the  slow  passage. 
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SIGNIFICANCE  AND  EXPLANATION 


Exothermic  chemistry  involving  Arrhenius  kinetics  is  of  practical 
importance  in  combustion  as  well  as  chemical  engineering.  Arrhenius  systems 
often  exhibit  multiple  steady  states;  a  typical  stationary  response  is  an 
S-curve.  From  such  a  steady  response  it  is  often  argued  that  the  system 
would  jump  from  a  weakly  reactive,  almost  extinguished  state  to  one  of 
vigorous  chemical  activity  as  a  control  parameter  (e.g.  the  Damkohler  number) 
is  increased  through  a  critical  value.  A  similar  jump  from  ignition  to 
extinction  is  implied  as  the  control  parameter  is  reduced  through  a  different 
critical  value.  Thus,  inherently  transient  pictures  are  inferred  from  steady 
solutions. 

The  purpose  of  this  study  is  to  provide  an  unsteady  description  of  the 
jump  phenomena.  Attention  is  focussed  on  situations  where  the  variation  of 
the  control  parameter  through  criticality  is  gradual.  Such  variations  might 
correspond,  in  practice,  to  a  slow  deterioration  in  the  activity  of  a 
catalyst  or  to  slow  increase  in  pressure.  The  asymptotic  analysis  also  makes 
use  of  large  activation  energy. 
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ARRHENIUS  SYSTEMS:  DYNAMICS  OF  JUMP  DUE  TO 
SLOW  PASSAGE  THROUGH  CRITICALITY 

* 

A.  K.  Kapila 

1.  INTRODUCTION 

Reactive  systems  involving  Arrhenius  chemistry  abound  in  several  areas 
of  practical  interest,  ranging  from  chemical-reactor  engineering  to  combus¬ 
tion.  Often,  such  systems  exhibit  multiple  steady  states.  Typically,  the 
equilibrium  plot  of  a  state  variable  y  (e.g.  the  maximum  temperature) 
against  a  control  variable  <j>,  measuring  the  reaction  rate,  is  an  S-curve 
(Figure  1) .  The  lower  branch  of  the  S  represents  the  low-conversion  or 
extinguished  state  while  the  upper  branch  corresponds  to  the  high-conversion 
or  ignited  state.  In  many  cases  both  of  these  branches  are  asymptotically 
stable  while  the  middle  branch  is  unstable,  the  exchange  of  stabilities 
occuring  precisely  at  the  turnaround  points  of  the  S.  The  steady  states 
and  their  stability  are  discussed  in  great  detail  by  Aris  [1,2]  in  the  con¬ 
text  of  diffusion  and  reaction  in  permeable  catalyst  pellets.  The  forth¬ 
coming  monograph  by  Buckmaster  and  Ludford  [3]  will  provide  several  relevant 
examples  from  combustion. 

This  paper  is  concerned  with  the  dynamic  response  of  the  system  to  slow 
variations  in  <|>.  Such  variations  may  be  due,  in  practice,  to  changes  in 
pressure  or  in  catalytic  efficiency.  Imagine,  for  example,  that  the  system 
is  initially  in  equilibrium  at  an  extinguished  state,  such  as  that  corre¬ 
sponding  to  point  A  in  Figure  1.  If  <j>  is  now  increased  gradually,  it  is 
expected  that  the  system  will,  more  or  less,  follow  the  lower  branch  until 
the  critical  point  (i.e.  the  ignition  point)  C  is  reached.  Then,  a  jump  to 
the  ignited  state  will  occur.  The  aim  of  this  paper  is  to  describe  this 
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transient  process.  A  similar  but  simpler  analysis,  not  presented  here,  de¬ 
scribes  the  jump  from  the  ignited  to  the  extinguished  state  due  to  slow 
passage  through  the  extinction  point  E. 

For  lumped  systems,  governed  by  nonlinear  ordinary  differential  equa¬ 
tions  of  type  dy/dt  =  F(y;4>),  the  response  due  to  slow  variations  in  <J> 
has  recently  been  examined  by  Haberman  [4].  He  looks,  in  particular,  at  a 
neighborhood  of  the  critical  point  C  and  shows  that  there  exists  a  transi¬ 
tion  region  which  connects  the  near-equilibrium  pre-critical  solution  (in 
which  4>  is  slowly  varying)  to  the  jump  solution  (in  which  <p  is  stationary 
at  $  ) .  For  those  cases  in  which  the  S-curve  is  parabolic  near  C,  Haberman 
shows  further  that  the  aforementioned  transition  region  is  governed  by  a 
Riccati  equation  whose  solution  explodes  in  a  finite  time.  We  shall  extend 
these  notions  to  show  that  an  entirely  analogous  situation  exists  for  the  dis¬ 
tributed  system,  governed  by  a  partial  differential  equation. 

This  paper  complements  earlier  investigations  by  Kapila  [5]  and  by  Kassoy 
and  Poland  [6,7]  on  the  dynamics  of  ignition  in  distributed  systems,  and  by 
Kassoy  [8]  on  similar  problems  in  lumped  systems.  There  the  emphasis  was  on 
evolution  at  fixed  <p ,  while  the  present  work  treats  time-varying  <p.  How¬ 
ever,  the  earlier  analyses  are  quite  relevant  here,  as  we  shall  see.  Slow- 
transition  problems  for  partial  differential  equations  have  also  been  con¬ 
sidered  by  Rubenfeld  [9] ,  but  only  for  bifurcation  points  rather  than  jump 
points.  Buckmaster's  treatment  of  slowly-varying  flames  [10]  also  falls  in 
the  former  category. 

Our  analysis  will  be  asymptotic,  depending  upon  the  largeness  of  two 
parameters,  one  characterizing  the  slow  variation  of  and  the  other  the 

activation  energy  of  the  reaction.  A  lumped  model  is  treated  first;  it  employs 
Haberman's  idea  in  a  slightly  amplified  form  and  lays  the  groundwork  for  the 
subsequent  analysis  of  a  distributed  model. 


2.  A  LUMPED  MODEL 


As  an  example  of  a  lumpted  system,  consider  the  Continuously  Stirred 
Tank  Reactor,  CSTR,  modelled  by  the  equation  [2] 

y^  =  A2(l+8-y)  exp(y-y/y)  +  1-y  ,  (2.1) 

where  y  is  the  temperature  and  t  the  time,  while  the  positive  parameters 
2 

A  ,6  and  y  are,  respectively,  the  Damkohler  number,  the  heat  of  reaction 
and  the  activation  energy.  It  is  convenient  to  define  a  reduced  Damkohler 
number  $  via  the  transformation 

A2  =  (e$y)  _1iJ)  .  (2.2) 

It  is  known  [2]  that  for  large  enough  y  the  equilibrium  plot  of  y 
against  $  has  an  S-shape  (Figure  1) ,  with  the  middle  branch  unstable  and 
the  extreme  branches  asymptotically  stable.  A  straightforward  analysis  of 
the  static  problem  (y^  =  0  in  (2.1))  in  the  limit  y  -+  °°  further  shows  that 
on  the  lower  branch  of  the  S  the  expansions 

y  =  1  +  y_1yi  +  0(y~2)  ,  <j>  =  ^  +  0(y_1)  (2.3) 

hold,  where  the  monotonically  increasing  function 

y1  =  G(((>1)  ,  0  <_  ^  <  1  (2.4) 

is  defined  implicitly  by 

y1  =  ^expfy  -1),  0  <  y1  _<  1  (2.5) 

and  is  shown  in  Figure  2.  The  ignition  point  C  corresponds  to 

=  y1  =  1  at  c  . 

Also,  for  future  use  we  note  that 


1 


(2.6) 


1/2  2  2 

G(ij>1)  =  1  -  [2(1-^1>]  7  +  -  (1-4>1)  +  0[(1-*1)  ]  as  $  - 

An  analogous  description  of  the  upper  branch  yields 

y=l+8-est  (2.7) 

D 

where  est  stands  for  terms  that  are  exponentially  small  as  y  °°,  and  the 
point  D  on  the  upper  branch  is  vertically  above  C. 

Of  interest  is  the  dynamic  behavior  of  y  as  (j  varies  slowly  through 
This  slow  variation  can  be  characterized  explicitly  by  introducing  a 
small  parameter  <5  and  a  slow  time  t,  where 

t  =  6t,  6  <<  1  , 

thereby  transforming  (2.1)  into 

6y  =  (e8y)  1<ji  (1+6-y)  exp(y-y/y)  +  1-y  .  (2.8) 

At  an  initial  time  t  =  <  0,  let  the  state  of  the  system  correspond  to 

point  A  in  Figure  1,  i.e.  let 

y(tA)  =  yA  ,  <MtA)  =  *A  - 

Then,  the  expansions  (2.3)  imply  that 

y(tA>  =  1  +  y  *y^  +  0(y  2)  /  <#>  ( t  A>  =  +  0(y  1)  .  (2.9) 

''  A  A 

Following  (2.3)  again,  we  let 

<f>  (t)  =  <t>1(t)  +  0  ( y”1)  (2.10) 

and  assume  that  <j>^(t)  is  a  smooth,  monotonically  increasing  function  which 
has  the  power  series  representation 

2 

$  (t)  =  1  +  t  +  0 (t  )  as  t  -*■  0  .  (2.11) 


4- 


This  specifies  t  =  0  to  be  the  instant  at  which  the  ignition  point  is 
reached;  there  is  no  loss  of  generality  in  setting  $|(0)  =  1.  The  goal  is 
to  obtain  an  asymptotic  solution  to  (2. 8-2. 9),  with  <j>^(t)  prescribed  above, 
in  the  limit  6  -*•  0,  y  -+  00 .  Several  different  regimes  need  to  be  distinguished 

2.1  Pre-critical  Solution 

Anticipating  that  the  solution  stays  close  to  unity  until  C  is  reached 
(i.e.  until  t  =  0) ,  we  substitute 

y  =  1  +  Y_1z(t)  (2.12) 

into  (2.8)  and  expand  for  y  -*■  “  to  get 

6zfc  =  <J>  ^  exp(z-l)  -  z  +  0  ( y  *)  ,  (2.13) 

where  (2.10)  has  also  been  employed.  Henceforth  the  analysis  will  concentrate 
on  the  limit 

6  >>  Y  1  .  (2.14) 

The  apparently  richer  limit  6  =  0(y  *)  can  also  be  treated,  at  the  simple 
expense  of  more  algebra,  but  was  found  to  yield  no  additional  significant 
effect . 

We  let  z  have  the  expansion 

z  =  z  +  6z2  +  o  ( <$)  (2.15) 

which,  when  fed  into  (2.13),  gives 

z1  =  4>1  exp(z  -1)  ,  (2.16) 

Z2  =  zi/(z1"1)  *  (2.17) 
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where  primes  stand  for  differentiation  with  respect  to  t.  Comparison  of 
(2.16)  with  (2.5)  shows  that  z^  has  the  same  form  as  the  static  solution 
(2.4),  i.e. 


z1  =  G ( <j>1  ( t)  ) 


(2.18) 


The  asymptotic  expressions  (2.6)  and  (2.11)  then  show  that,  as  expected, 

(2.18)  is  valid  only  for  t  <  0.  In  fact,  we  find  that 

1/2  2  1/2 

z  1  =  1  ~  (-2t)  '  -  -  t  +  0(-t)  '  as  t  -*•  0  .  (2.19) 

Once  z^  is  known,  (2.17)  determines  z^  and  we  get  the  asymptotic  behavior 

z2  =  ( 2t)  ” 1  +  j(-2t)_1/2  +  0(1)  as  t  -<•  0~  .  (2.20) 

The  above  expression  exhibits  even  more  dramatically  the  breakdown  of  the 
pre-critical  expansion  (2.15). 

A  word  about  the  initial  conditions  (2.9)  is  in  order.  These  imply  that 

z^(tft)  =  y^  and  z2(tfl)  =  °*  While  z  ,  as  given  by  (2.18),  satisfies 
A 

this  requirement,  z^,  as  given  by  (2.17),  does  not,  in  general.  This  dis¬ 
crepancy  is  easily  remedied  by  means  of  a  thin  initial  layer  at  tft  in  which 
the  relevant  time  is  the  fast  time  x.  In  this  layer,  the  effect  of  the 
initial  condition  on  z^  decays  exponentially  rapidly.  In  fact,  even  initial 
conditions  off  the  lower  branch  can  be  accomodated  in  this  manner  provided, 
of  course,  that  the  initial  point  lies  within  the  domain  of  attraction  of  the 
(asymptotically  stable)  lower  branch. 


2.2  Transition  Solution 

Further  development  of  the  solution  occurs  on  a  new  time  scale  s,  de¬ 
fined  by  the  stretching 


-6- 


(2.21) 


1/3.2/3 

t  =  -2  6  s 


-1/3 

The  new  time  is  short  on  the  t  scale  but  is  still  long,  0(6  ),  on  the 

t  scale.  Equations  (2.11)  and  (2.13)  transform  into 


4>1  =  1  -  21/362/3s  +  0(64/3) 


(2.22) 


and 


-2_1/361/3z  =  [1-21/362/3s  +  0(64/3)  +  0 ( y"1 ) ] exp (z-1)  -  z  +  0(y_1)  ,(2.23) 


and  the  new  expansion  for  z  is  taken  to  be 


z  =  1  +  (s)  +  62//3v2(s)  +  0(6) 


(2.24) 


It  is  clearly  shown  by  (2.22)  and  (2.24)  that  the  transition  solution  is  re¬ 
stricted  to  a  small  neighborhood  of  the  critical  point  C.  Substitution  of 
(2.24)  into  (2.23)  yields  the  following  Riccati  equation  for  v  : 


v1=21/3(22/3s-iv2)  ,  s< 


(2.25a) 


The  condition 


2/3  ,1/2  1  , 

v,  -2  (s  +  -— )  as  s 
i. 


(2.25b) 


comes  from  matching  (2.24)  with  the  pre-critical  solution.  The  solution  to 
(2.25a,b)  is  (see  [4] ) 


v1  =  22/3Ai’ (s)/Ai(s)  , 


(2.26) 


where  Ai  is  the  Airy  function.  This  solution  is  valid  only  for  s  >  s  , 


where  sQ  =  -2.3318  is  the  first  zero  of  Ai(s),  and  there  it  has  the 
explosive  behavior 


2/3  -11  + 

V1  =  2  Ks-s0)  +  3  s0(s-s0)]  as  s  -  sQ 


(2.27) 
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This  singularity  immediately  signals  the  breakdown  of  the  transition  solution 
and  suggests  that  the  solution  is  attempting  to  deviate  from  criticality 
farther  than  is  allowed  by  the  expansion  (2.24).  Before  proceeding  to  the 
next  regime,  however,  it  is  desirable  to  compute  v  as  well.  The  problem 
for  v^  is  found  to  be 


whence 


v^  =  "2l^3(v1v2  +  i  vi  "  2’L//3sv1)<  s  <  °° 
21,/3  1  -1/2 

v2  ->  —  (2s  +  —  s  )  as  s  -*■«*>  , 


(2.28) 


V2  "  2 


1/3 


s  +  {Ai ( s) }  2  f  {Ai(x)}2dx 
-  s 

00  ' 
+  y{Ai(s) }  2  /  {  Ai ' (x) }3{  Ai  (x)  }  3dx 


(2.29) 


and  it  has  the  asymptotic  behavior 

24/3  -2  -2  + 
v2 - J  ^S~S0)  ?'n<s*SQ)  +  A0(s-s0)  ,  as  s  -  sQ 


(2.30) 


The  constant  Aq  is  given  by 


AQ  =  21/3{Ai' (sQ) }~2 


|4  {  J  {Ai'  (x)  }3{Ai(x)  }  1dx 


0 


+  j  ({Ai'(x)}  {Ai (x) }  -  {Ai' (sQ) }2(x-sQ)  ) dx 

so 


+  {Ai ' (sQ) }Z£n (-sQ)  \  +  J  Ai2 (x) dx 


(2.31) 
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2.3  Post-Critical  Solution 

The  manner  in  which  (2.24)  breaks  down  suggests  the  stretching 

s  =  sQ  -  2_1/361/3(p-p0)  (2.32) 

where  the  new  time  scale  p  is  of  the  same  order  as  the  fast  time  x,  and 
the  shift  Pq,  assumed  to  be  o(6  ^/'^),  is  to  be  determined.  The  expression 
v2.22)  reduces  to 

=  1  -  21/362/3sq  +  6(p-p0)  +  0(64/3)  , 

indicating  that  <p  remains  close  to  (but  larger  than)  the  critical  value, 
while  (2.23)  transforms  into 

2/3  -1 

z  =  [1  +  0(6  )]exp(z-l)  -  z  +  0(6)  +  0(y  ) 

P 

We  note  that  the  unsteady  term  z^  now  appears  to  leading  order.  Allowing 
z  to  have  the  expansion 

z  =  w  (p)  +  o(l)  ,  (2.33) 

w^  is  found  to  satisfy 

w|  =  exp (w^-1 )  -  w  ,  p  >  -»  , 

yielding  the  implicit  solution 

Wi 

/  (exp(x-l)  -  x]  dx  =  p-p  .  (2.34) 

2  u 

The  integration  constant  pQ  is  to  be  determined.  It  can  be  shown  that 

,24  A1 

W1  ~  1  ~  p"  ~  — 2  +  — 2  as  p  , 

3p  p 

where 


-9- 


+  2 


*1  "  I  in  2  -  4  -  2>'0 


J  l  (exp(x-l) -x}  J‘-  2  (x-1)  2  +  -(x-1)  ^)dx  .  (2.35) 


Matching  with  the  transition  solution  (2.24)  yields 


»0  =  -9Jn(5)  ' 


confirming  the  expectation  that  =  o  (6  and 


\  =  22/3A0  +  |  ^  2  . 


With  known  from  (2.31),  substitution  of  the  above  expression  into  (2.35) 

will  determine  p  . 

At  this  stage  it  is  convenient  to  revert  to  the  original  dependent 
variable  y.  Substitution  of  (2.33)  into  (2.12)  yields  the  expansion 


y  =  l  +  y  1w1(p)  +  o(y  1)  / 


(2.36) 


indicating  that  so  far  y  has  stayed  close  to  unity.  However,  larger  depart¬ 
ures  from  unity  are  imminent,  because  the  solution  (2.34)  for  explodes 

in  a  finite  time.  In  fact,  (2.34)  shows  that  increases  monotonically 

and  has  the  behavior 


w,  =  -5,n(p  -  p)  +  1  +  o(l)  as  p  •*  p 

]  00  0 


p  =  Pn  +  f  [exp(x-l)  -  x]  dx 

oo  0 


(2.37) 


(2.38) 


2.4  Jump  Solution 

The  singularity  (2.37)  implies  that  (2.36)  is  no  longer  valid,  and  that 
further  development  of  the  solution  takes  place  on  the  exponentially  fast  time 
scale  ;,  defined  by  the  nonlinear  stretching 


p 


p  -  exp(-vo) 

CO 


a  >  0 


On  the  o  scale,  (2.10)  reduces  to 


=  1 


„l/3 ,2/3 
2  6  s, 


+  <5  (p_  -  p. 


e'Y°] 


0(64/3) 


0  ( Y_1) 


indicating  that  $  undergoes  only  exponentially  small  variations  from  a 

2/3 

fixed  value  in  a  6  -neighborhood  of  criticality.  Therefore  Kassoy's  super¬ 
critical  analysis  [8]  ,  carried  out  for  fixed  <)>,  becomes  relevant.  We  shall 
not  repeat  the  details  here;  suffice  it  to  say  that  the  solution  has  the 
expansion 


y  =  (l-o)  1  +  0 ( y  3) 

and  that  y  reaches  the  upper-branch  value  1+6  (cf.  (2.7))  at  o  =  B/(l+6) 
i.e.  at 


t  =  -21/362/3sn  +  6  Ip  -  p  -  exp{-YB/(l  +  B) })  - 

0  or  (J 

We  note  that  during  this  exponentially  rapid  ascent,  the  governing  parameter 
is  y<‘  hardly  any  role  is  played  by  6.  However,  once  the  upper  branch  has 
been  reached,  further  evolution  of  the  solution,  involving  slow  variation 
about  the  upper-branch  equilibrium,  will  be  governed  by  a  6-expansion 
analogous  to  (2.15). 
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3.  A  DISTRIBUTED  MODEL 


As  an  example  of  distributed  systems,  consider  the  porous  catalyst 
pellet,  for  unit  Lewis  number  and  in  the  symmetric  slab  geometry,  governed  by 
the  equations  [1] 

2 

y^  =  y  +  X  (1+6-y) exp(Y-y/y} ,  0  <  x  <  1  ,  (3.1a) 

y^(0,T)  =  0,  y(l,T)  =  1  ,  (3.1b) 

with  appropriate  initial  condition.  The  symbols  have  the  same  meanings  as  in 
Section  2;  the  new  independent  variable  x  denotes  the  spatial  coordinate. 

It  is  again  convenient  to  eliminate  X  in  favor  of  <j>  through  the  relation 

X2  =  (0y)~14>  •  (3.2) 

The  relevant  static  problem,  which  is  a  classic  in  the  chemical  engineer¬ 
ing  literature  [1] ,  was  recently  studied  by  Kapila  and  Matkowsky  (see  the 
Appendix  of  [11]  )  in  the  limit  y  •*  <=°.  Schematically  the  response  diagram, 
now  a  plot  of  y  at  x  =  0  against  4>,  is  the  same  as  in  Figure  1,  with  the 
same  stability  properties  [2].  On  the  lower  branch  the  expansions  (2.3)  still 
prevail,  leading  to  the  problem 

y1  +  ^  exp  y^  =  0,  y^  (0)  =  y^O)  =  0  (3.3) 

XX  X 

whose  solution 

yx(x)  =  H  ( x ;  x )  (3.4) 


can  be  represented  parametrically  by  the  equations 


H  =  2  in  (cosh  a  sech(ax)], 


2  2 

2a  sech  a,  0  <  a  <  a 


(3.5) 


(Schematically,  Figure  2  also  represents  the  graph  of  y^tO)  against  $  . ) 
The  critical  value  ac  is  given  by 
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I 


a  tanh  a 
c  c 


1 


(a  ~  1.2) 
c 

whence 

<b1  *>  0.88,  y  (0)  *  1.187  . 

c  c 

The  analysis  in  [11]  also  shows  that  at  the  point  D  on  the  upper  branch 
(Figure  1) , 

yD<0)  =  1+0  -  est  .  (3.6) 

In  fact,  y  (x)  =  1+8  -est  throughout  the  domain,  except  for  an  exponentially 
thin  layer  at  x  =  1. 

We  shall  find  that  the  similarity  between  the  static  behaviors  of  the 
lumped  and  the  distributed  systems  also  extends  to  the  dynamic  responses  due 
to  slow  variation  of  <j>.  Before  pursuing  that  question,  however,  we  digress 
for  a  brief  look  at  the  linear  stability  problem  for  the  lower  branch.  This 
problem  can  be  shown  to  yield  the  following  leading-order  eigenvalue  problem 
in  the  limit  y  + 

U"  +  ($  exp  y  )U  -  AU  =  0,  U' (0)  =  U ( 1)  =  0  .  (3.7) 


(Equation  (3.7)  can  be  derived  most  simply  by  reinstating  the  time-derivative 
y^  in  (3.3),  linearizing  about  the  steady  state,  and  seeking  solutions  of 

T 

the  type  U(x)  exp(Ax).)  Clearly,  the  eigenvalues  and  the  corresponding 

eigenfunctions  lh(x)  of  this  Sturm-Liouville  problem  depend  upon  <f>  .  For 
later  use,  we  note  two  properties  of  (3.7).  First,  all  A^  are  negative 
for  *  4^  >  as  a  consequence  of  the  asymptotic  stability  of  the  lower 


branch.  In  particular,  therefore,  A  =  0  is  not  an  eigenvalue  for  $  <  <(> 

c 

Second,  A  =  0  is.  an  eigenvalue  at  the  critical  point  <Ji  =  <J>  ,  because 

c 

that  is  where  the  exchange  of  stabilities  occurs.  The  corresponding  eigen¬ 
function  is  given  by 
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U  (x)  =  1-a  x  tanh  (a  x) 
20  c  c 


(3.8a) 


1 


Incidently,  another  (linearly  independent)  solution  of  the  differential  equa¬ 
tion  in  (3.7)  is 

uin(x)  =  —  tanh  (a  x)  .  (3.8b) 

10  a  c 

c 

Turning  now  to  the  dynamic  problem,  we  again  employ  the  slow  time 
t  =  St,  let  y  =  y(x,t)  and  rewrite  (3.1a,b)  as 

<5yt  =  yxx  +  (By)  1<P  exp(y-y/y)  ,  0  <  x  <  1  ,  (3.9a) 

y  <0,t)  =  y (1 ,t)  =  1  ,  (3.9b) 

and  again  take  the  initial  condition,  at  t  =  t  <  0,  to  be  that  correspond¬ 
ing  to  point  A  on  the  lower  branch  (Figure  1).  The  prescription  (2.10)  for 
still  holds,  but  now  we  let 

<|>  (t)  =  <p1  [1  +  t  +  0 ( t2 )  ]  as  t  ->  0  .  (3.10) 

c 

The  corresponding  behavior  of  the  parameter  a,  appearing  in  (3.5),  is 

a  (t)  =  a  -  (-t)1/2  -  ~  t  +  0  [  (-t) 3/2  ]  as  t  +  O-  .  (3.11) 

c  3a 

c 

This  result  will  be  found  useful  in  future  computations. 

As  in  the  lumped  example,  the  evolution  of  the  solution  is  again  followed 
through  several  stages. 


3.1  Pre-critical  Solution 

Expecting  y  to  stay  close  to  unity  prior  to  criticality,  we  substitute 


y 


1  + 


-1 

Y 


z (x, t) 


(3.12) 


into  ( 3 . 9a ,b)  to  get 
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(3.13) 


5zt  = 


0(Y_1)  , 


z  (0,t) 

x 


Z(l,t)  =  0 


For  6  >>  y  ,  we  assume  the  expansion 

z  =  z^  +  6z^  +  o(6)  (3.14) 

and  thereby  obtain,  for  z^,  the  problem 


z^  +  4>  exp  z  =  0,  z1  ( 0 , t)  =  z  (l,t)  =  0  . 

XX  X 

Comparison  with  (3.3)  yields  the  pseudo-steady  solution 

z  L  (x ,  t )  =  H(x;<)>  (t))  .  (3.15) 

Application  of  the  asymptotic  expressions  (3.10)  and  (3.11)  into  the  definition 
(3.5)  of  H  then  shows  that 


2  1/2 

z  (x)  -  —  (-t)  u  (x)  +  Oft) 
l  a  2u 

c  c 


as 


t  -*•  0 


(3.16) 


where 


z,  (x)  =  2  £n[cosh  a  sech(a  x) ]  , 

1  c  c 

c 


(3.17) 


and  u2Q(x)  has  already  been  introduced  in  (3.8a).  The  behavior  (3.16), 
analogous  to  the  lumped  result  (2.19),  demonstrates  that  (3.15)  is  not  valid 
beyond  t  =  0. 

It  is  found  that  z satisfies 


z2  +  (<j>^  exp  z^)z2  ~  zi  ’  z2  <0.t)  =  z2(l,t)  =  0,  t  <  0  .  (3.18) 

XX  t  X 

This  nonhomogeneous  problem  yields  a  unique  z2,  as  shown  by  the  following 

argument.  The  homogeneous  problem  corresponding  to  (3.18)  is  identical  to 

the  eigenvalue  problem  (3.7)  for  A  =  0  and  for  <(>  <{(>  ,  and  therefore  has 

c 


I 
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only  the  trivial  solution  (see  the  remark  about  zero  not  being  an  eigenvalue 


for  $  <  $  ,  made  following  (3.7)).  Hence  z^  is  unique,  and  has  the 

c 

representation 


z  ^x,t)  =  a(t)u  (x,t)  +  /  [u  (x,t)u  (C.t) 
2  z  o  1  2 


u1(C,t)u2(x,t) ]z1  (£,t)d£  . 


:3.19) 


Here, 

,  1 

a ( t)  =  -[u  (l,t)]  /  [u  (l,t)u  (£,t)  -  u. (e,t)u_(l,t)]z.  (C.t)dC 

0  1  2  12  1t 

and  the  u^(x,t)  satisfy  the  homogeneous  differential  equation 

Ui  +  exp{z^ (x,t) }]u^  =  0  (3.20a) 

XX  X  J 

and  the  initial  conditions 


ux  (0,t)  =  u2(0,t)  =  lj  vi^Opt)  =  u2  (0,t)  =  0  .  (3.20b) 

X  X 

(The  remarks  made  at  the  end  of  Section  2.1  about  satisfaction  of  initial 
conditions  by  z^  and  z2  apply  here  as  well.)  The  asymptotic  nature  of 
u^  and  u2  is  investigated  in  the  Appendix  and  the  results  are 

1/2 

Uj.tXjt)  =  u1Q(x)  +  (-t)  u  (x)  +  °(t)  »  (3.21a) 

1/2 

u2(x,t)  =  u2Q(x)  +  (-t)  u2i(x)  +  as  t  -+  0  .  (3.21b) 

The  functions  u^  and  u^  have  been  defined  in  the  Appendix,  while  u-^0* 
u2Q  were  introduced  in  (3.8a,b) .  With  these  expressions  in  hand,  (3.19)  can 
be  expanded  for  small  t  and  the  result  is 

z2(x,t)  =  — ■  (t)  1u2Q(x)  +  0[(-t)  1/2]  as  t  -*  0~  ,  (3.22) 

3a 

c 


/ 
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which  is  analogous  to  (2.20).  Thus,  the  pre-critical  expansion  (3.14)  becomes 
disordered  in  much  the  same  way  as  its  lumped  counterpart  (2.15)  did,  and  we 
are  led  to  the  transition  solution. 

3.2  Transition  Solution 
The  transition  time  scale  s  is  now  defined  by 

t  =  -62/3b-1s  ,  (3.23) 

where  the  0(1)  constant  b  will  be  chosen  later,  in  a  way  that  simplifies 
algebra.  When  the  above  stretching,  and  the  assumed  transition  expansion 

z  =  z^  (x)  +  61//3v^(x,s)  +  62//3v2(x,s)  +  6v^(x,s)  +  o(6)  ,  (3.24) 

c 

are  fed  into  (3.13),  the  following  equations  for  the  emerge: 


V1 

+ 

<♦1 

exp 

z!  )V1  =  0  ' 

(3.25a) 

XX 

c 

c 

-1  1  2. 

v2 

+ 

(*1 

exp 

2!  )v2  =  -bvi  " 

(♦l 

exp  z1  )  v-b  s  +  2"  vi^  ' 

(3.25b) 

XX 

c 

c  s 

c 

c 

V3 

+ 

(♦l 

exp 

zi  )V3  =  -bV2  - 

(♦l 

,  ,  13  -1 

exp  z 1  ) (v  v2  +  g  vx  "  b  svx) 

.  (3.25c) 

XX 

c 

c  s 

c 

c 

Each 

V  . 

1 

is. 

of  < 

course,  subject 

to  the  boundary  conditions 

v^  (0,s)  =  v^ 

(l,s) 

=  0;  i  =  1,2,3  . 

(3.26) 

x 

The  problem  for  can  be  solved  immediately  to  yield 

Vl(x,s)  =  f1(s)u2Q(x)  (3.27) 


where  the  "amplitude  function"  f^s)  is  to  be  determined  by  requiring  that 
the  nonhomogeneous  equation  (3.25b),  subject  to  (3.26),  have  a  solution.  This 
requirement  leads  to  the  orthogonality  condition 
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1 


[-bf^(s)  u90(x)  -  exp  ( x )  }  { -b  1s  +  f  ^  (s )  u^q  (x)  }  ]  u2Q  (x)  dx  , 

c  c 


which,  upon  t  he  evaluation  of  appropriate  integrals,  reduces  to  the  differential 
equation 


bf'(s)  =  3b  1s  -  ja2f^(s),  s  < 
1  4  c  1 


The  condition 


f,  - - 


2  -1/2  1/2  b  -1 

~  b  s  -  — -  s)  as  s  -*■ 
a  2 

c  3a 

c 


is  provided  by  the  pre-critical  solution  upon  matching.  Thus,  the  familiar 
Riccati  problem  emerges  again  (cf.  (2.25a,b)).  In  fact,  the  choice 

b  =  (9a2/4)1/3 
c 

leads  to  the  solution 

2  -2/3 

f x ( s)  =  3(9ac/4)  Ai'(s)/Ai(s)  ,  (3.28) 

thereby  determining  v1  completely,  once  (3.28)  is  substituted  into  (3.27). 

It  is  now  possible  to  write  down  the  solution  to  (3.25b)  and  (3.26).  The  re¬ 
sulting  expression  for  v 2  will  contain  f2(s)u  Q(x)  as  the  complementary 
function.  By  appealing  to  the  solvability  requirement  on  the  v^-equation, 
(3.25c),  and  by  matching  with  the  pre-critical  solution,  we  can  show  that  f^ 
satisfies  a  differential  problem  (analogous  to  that  in  (2.28))  which  can  be 
solved.  The  relevant  calculations  were  made,  because  they  are  needed  to 
determine  the  leading-order  solution  at  the  next  stage.  However,  although 
straightforward,  they  are  too  cumbersome  to  report  here. 

As  s  approaches  s^ ,  f^  explodes  (see  (2.27)),  leading  us  to  the 
next  stage  of  solution  development. 
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3.3  Post-critical  and  Jump  Solutions 
Following  the  treatment  of  Section  2.3,  the  new  time  variable  p  is 
defined  by 

1/3 

s  =  sQ  -  b6  (p-p  ) 

or  equivalently  by 

t  =  -62/3b_1s0  +  6(P-P0)  , 

where  pQ  is  the  yet  unknown  shift  as  before,  and  z  is  assumed  to  have  the 
expansion 

z  =  w1(x,p)  +  o(l)  -  (3.29) 

Then,  (3.13)  provides  the  leading-order  problem 

w^  =  w^  +  <fi^  exp  w^,  0  <  x  <  1 ,  p  >  -°°  ,  (3.30a) 

xx  c 

w  (0,p)  =  w1(l,p)  =  0  ,  (3.30b) 

x 

for  which  the  initial  condition 

W1  =  Z1  *x* - ~ U20^X^  +  as  P  "*•  ~°°  (3.30c) 

c  3a  p 

c 

comes  from  matching  with  two-terms  of  the  transition  expansion  (3.24).  In 
order  to  determine  the  shift  pQ  and  to  fix  the  origin  of  p  in  the 
problem  (i.e.  to  fix  the  counterpart  of  the  constant  pQ  appearing  in  (2.34)), 
matching  with  three  terms  of  (3.24)  is  needed,  and  this  is  where  the  knowledge 
of  v2(x,s)  is  required. 

We  now  stress  two  aspects  of  w  .  First,  it  evolves  at  a  constant  value 

ip  of  Second,  it  measures  the  departure  of  y  from  unity  on  the 

C-1 

0(y  )  scale;  this  is  shown  by  the  asymptotic  expression 
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(3.31) 


y  =  1  +  Y  1«1 (x,p)  +  o ( y  1) 

obtained  by  substituting  (3.29)  into  (3.12).  Therefore,  w^  plays  the  same- 
role  as  the  induction-period  solution  of  the  author's  earlier  paper  [5],  where 
transition  from  an  extinguished  to  the  ignited  state  is  studied  at  fixed  1 ; 
also  see  [6).  In  fact,  as  we  shall  see,  the  entire  analysis  of  [5]  applies 
here,  as  did  Kassov's  fixed  tp  analysis  [8]  in  the  lumped  case. 

Figure  3  shows  the  numerical  solution  of  (3.30a,b,c),  obtained  by 
using  the  package  PDECOL  [12]  .  This  package  employs  B-splines  for  spatial 
discretization,  and  then  integrates  the  resulting  ordinary  differential  equa¬ 
tions  with  a  Gear  solver.  Integration  was  begun  at  a  large  negative  value 
of  p.  We  find  that  initially  the  solution  develops  slowly,  but  then  w^ 
begins  to  rise  rapidly  near  x  =  0  while  variations  continue  to  be  leisurely 
in  the  rest  of  the  domain.  Eventually,  at  a  definite  p  =  p  ,  w  (0,p) 

OO 

becomes  unbounded  (compare  with  Figure  3  of  [5]). 

The  explosive  behavior  of  w^  in  the  thin  boundary  layer  at  x  =  0 
can  be  analyzed;  this  was  done  in  [5]  and  a  logarithmic  singularity  (cf.  (2.37)) 
was  revealed.  This  singularity  marks  the  breakdown  of  (3.31),  and  further 
development  is  exponentially  fast,  as  it  was  in  the  lumped  case.  For  full 
details  the  reader  is  referred  to  [5] .  Summarizing  briefly,  a  rapidly  shrink¬ 
ing  hot  spot  of  growing  intensity  developed  at  x  =  0.  When  y  reaches  the 
valued  1+6  there,  the  hot  spot  detaches  from  the  left  boundary  and  moves 
into  the  interior  of  the  domain  as  a  well-defined  reaction  wave,  still  travel¬ 
ling  exponentially  fast  and  leaving  behind  a  burnt  zone  at  constant  y  =  1+3. 
Just  before  the  right  boundary  is  reached,  the  wave  quickly  comes  to  rest  (on 
the  exponentially  fast  scale)  to  accomodate  the  boundary  condition  at  x  =  1. 

In  other  words,  the  jump  to  point  D  in  Figure  1  is  now  complete  (see  (3.6)). 
Further  movement  along  the  upper  branch  will  again  be  governed  by  the  slow 
variable  t,  much  in  the  manner  of  the  pre-critical  solution. 
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4 .  CONCLUDING  REMARKS 


The  asymptotic  analysis  has  concentrated  on  the  behavior  of  the  systems 
as  the  parameter  <j>,  a  measure  of  the  rate  of  reaction,  passes  slowly  through 
the  ignition  criticality.  It  is  shown  that  the  solution  follows  the  lower 
branch  at  the  slow  scale  6t  (where  t  is  the  reference  scaie,  i.e.  the 
residence  time  for  lumped  systems  and  diffusion  time  for  distributed  systems) , 
passes  through  criticality  at  the  slightly  faster  scale  and  goes 

through  the  initial  stage  of  the  jump  at  the  scale  t.  Most  of  the  jump,  how¬ 
ever  occurs  at  an  exponentially  fast  scale  (in  y) ,  and  the  subsequent  travel 
along  the  upper  branch  occurs  again  at  the  slow  scale  6x.  The  similarities 
between  the  responses  of  the  lumped  and  distributed  systems  are  clearly 
demonstrated. 

From  a  practical  viewpoint,  the  ignition  problem  considered  here  is 
probably  of  greater  relevance  in  combustion,  where  increase  of  if>  may  be 
associated  with  increase  of  pressure.  On  the  other  hand,  in  chemical 
engineering  the  extinction  problem  might  have  more  significance,  where 
decrease  in  41  could  be  the  result  of  a  decline  in  the  efficiency  of  a  catalyst. 
A  similar  treatment  would  apply  there,  but  the  extinction  jump  would  be  much 
less  dramatic. 
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APPENDIX 


Here  we  are  concerned  with  the  solutions  u.(x,t),  i  =  1,2,  of  the  euua- 

1 


u  +  F(x,t)  =  0 

XX 


subject  to  the  initial  conditions 


u1  (0  ,t)  =  1,  (0 ,t)  =  0;  u2  (0, t)  =  0,  u2(0,t)  =  1 

X  X 


The  function  F  appearing  in  (A.l)  is  defined  by 


F (x , t)  =  ij>  (t)  exp[z^  (x,t)  ] 


which,  from  (3.15)  and  (3.5),  is  seen  to  be 


F(x,t)  =  2a^  sech^(ax)  .  (A. 3 

Of  particular  interest  is  the  behavior  of  u.^  and  u2  as  t  -*  0  .  Firs 

2  2 

we  note  that  at  t  =  0,  a  =  a  so  that  F(x,0)  =  2a  sech  (a  x) .  Then,  it  is 

c  c  c 

easily  checked  that  (A.l)  and  (A. 2)  have  the  solutions 


u  (x,0)  =  u  (x)  =  —  tanh  (a  x)  , 
1  10  a  c 

c 


u_(x,0)  =  u.^x)  =  1  -  a  x  tanh  (a  x) 
2  20  c  c 


These  functions  have  already  appeared  in  (3.8a,b).  Now,  in  view  of  the 
asymptotic  expression  (3.11)  for  a,  (A. 3)  gives 

2  2  2  1/2 

F(x,t)  =  2a  sech  (a  x)  [1  -  —  (-t)  u„„(x)  +  0 ( t )  ]  as  t  -+  0  ,  (A. 4 

c  c  a  20 

c 


suggesting  the  expansions 


u.  (x,t)  =  u.~(x)  +  (-t)  u.,(x)  +  0(t)  as  t  -*•  o  ;  i  =  1,2 


When  (A. 4)  and  (A. 5)  are  substituted  into  (A.l)  and  (A. 2),  we  find  that  the 


ir^(x)  satisfy 


2  2  2 
u"  +  2a  sech  (a  x)u.,  =  4u.  u_„a  sech  (ax);  l 
ll  c  c  ll  lO  20  c  c 


=  1,2  , 


u!  (0)  =  u.  .  (0)  =  0  , 

ll  ll 


leading  to  the  solutions 


u  (x)  =  /  4a  u  (£)u  (C)sech2(a  ?)[u  (x)u  (£)  -  u  (£) u  (x) ]d£;  i  =  1, 
ll  •_  c  lO  20  c  10  20  10  20 


Upon  evaluation  of  the  appropriate  integrals,  we  find  that 


1  2 
u  (x)  =  — —  [tanh  (ax)  -ax  sech  (ax)]  , 

11  2  c  c  c 

a 

c 


u  , (x)  =  x[tanh  (ax)  tax  sech  (ax)] 
21  c  c  c 


FIGURE  CAPTIONS 


Figure  1.  The  steady-state  response 
Figure  2.  The  lower  branch 


Figure  3.  The  numerically-computed  plot  of  w^(x,p) 
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